Subordination model of anomalous diffusion leading to the two- 
power-law relaxation responses 
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Abstract. - We derive a general pattern of the nonexponential, two-power-law relaxation from 
the compound subordination theory of random processes applied to anomalous diffusion. The 
subordination approach is based on a coupling between the very large jumps in physical and 
operational times. It allows one to govern a scaling for small and large times independently. 
Here we obtain explicitly the relaxation function, the kinetic equation and the susceptibility 
expression applicable to the range of experimentally observed power-law exponents which cannot 
be interpreted by means of the commonly known Havriliak-Negami fitting function. We present 
a novel two-power relaxation law for this range in a convenient frequency-domain form and show 
its relationship to the Havriliak-Negami one. 



Introduction. — Many studies have been reported on 
the phenomenon of nonexponential, power-law relaxation 
which is typically observed in complex systems such as 
dielectrics, ferroelectrics, polymers and so on (see [TJ|2 
and references therein). The main feature of such sys- 
tems is a strong (in general, random) interaction between 
their components in the passage to a state of equilibrium. 
Therefore, theoretical construction of an "averaged" ob- 
ject representing the entire relaxing system is not a simple 
problem; it requires application of advanced probabilistic 
tools. One of them uses a randomization in the parame- 
ters of distributions that describes the relaxation rates in 
disordered systems. With regard to the dielectric relax- 
ation, each individual dipole in a dielectric system relaxes 
exponentially, but their relaxation rates are different and 
obey a probability distribution (continuous function) [3|. 
This method is successive for getting many empirical re- 
sponse laws and their classification, but it sometimes be- 
comes enough complicated to interpret their interrelations 
and to derive kinetic/diffusion equations. An alternative 
approach, applied to analysis of diffusion processes under- 
lying the macroscopic dynamics, is based on subordina- 
tion of random processes (see, for example, [4] and refer- 
ences therein). The anomalous diffusion process y[i7(t)] 
is obtained by the time clock randomization of a random 
(parent) process Y{t) by means of another random pro- 



cess U{t) called the directing process. The latter process 
is also often referred to as the randomized or operational 
time [5]. As it has been shown in [6,7 , the subordina- 
tion method is also useful for analysis of the relaxation 
processes. However, up to now the approach had some 
difficulties in interpretation of the Havriliak-Negami (HN) 
law, one of the most general description of relaxation data. 
Only recently the effective picture of the HN relaxation, 
based on subordination, has been found in [8 . It requires 
introduction of compound subordinators to the study of 
the anomalous diffusion. The compound subordination 
procedure relates one subordinator to another in such a 
way that the temporal evolution of the parent random 
process has a different scaling on short and long times. 

Coupling between the a-stable process as the 
real time and its inverse as an operational time. 

— The inverse a-stable process is actually the left-inverse 
process of the a-stable process (see Fig. [T]). This means 
that the a-stable process T(r) and its inverse St = inf {r > 
I T{t) > t} fulfill the relation St(t) = t, while T{St) > t 
holds. The random variable St corresponds to the first 
passage time of the strictly increasing process T{t) above 
t. To find the probability density distribution (pdf) of 
the compound process T{St), we start with a sum of inde- 
pendent and identically distributed heavy-tailed random 
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Fig. 1: (Color online) Trajectory of T{t) and St with a = 0.75. 

variables T^, namely J7„ = X^ILo with Tq = 0. The pro- 
cess Nf — max{n G N | [/„ < i} is inverse to J7„. From 
this definition it follows directly the inequality 



UNt <t<U 



Nt + l 



for t > . 



(1) 



According to the Feller's book [5], the renewal theory tells 
us how the pdfs of U and Un^^i behave in the long time 
limit. The random variable Uj^^/t tends (in distribution) 
to Y with the pdf equal 



sin(7ra) ^ 



{l-xY 



< a; < 1 , 



and UNt+i/t tends to Z which pdf reads 
sin(7ra) 



(2) 



(3) 



The functions (x) and (x) correspond to special cases 
of the well-known beta density [5^. We note that the pdf 
oip^ (x) concentrates near and 1, whereas p^ (x) concen- 
trates near 1, where both tend to infinity. The moments 
of the random value Y are finite and can be calculated 
directly from the pdf However, even the first moment 
of Z is infinite. 

A passage from the discrete process Tj to the continuous 
limit T(r) allows one to rewrite the inequality ([T]) into a 
form 

T-{St) <t<T{St) for t>0, (4) 

underestimating or overestimating the real time t. In this 
case we obtain an analogous to and ([3]) result. In the 
long hmit the pdfs of T-{St) and T{St) read 



p+{t,z) 



sm TTa 
sin TTa 



y''-\t~y)-^, 0<y<t, (5) 



y-ir(y-t)-", y>t. 



(6) 



Their moments are computed directly from the moments 
of Y and Z due to the relations T-{St) = tY and T{St) = 



tZ, where ~ denotes the equality in distribution. Thus, 
the process T~{St) has finite moments of any order, while 
T(S't) gives us even no finite the first moment. It should 
be noticed that the processes T"(S't) and T{St) evolve to 
infinity with growing t. They show a coupling between the 
very large jumps in physical and operational times. 

Universality of relaxation in two-state systems. — 

The simplest interpretation of relaxation processes uses 
the concept of a system of independent exponentially re- 
laxing objects (for example, dipoles) with different (inde- 
pendent) relaxation rates [9]. The systems, following this 
law (called Debye's), can easily be described as a two- 
state system. Let N be the common number of dipoles in 
a dielectric system. If is the number of dipoles in the 
state t, is the number of dipoles in the state ], so that 
N — + N^. Assume that ioi t = the system is stated 
in order so that the states t dominate, namely 



N^(t ^ 0) 
N 



= %(0) = 1, 



m ^ 0) 

N 



= "a(0) = 0. 



where is the ratio of dipoles in the state f and ni in the 
state Denote the transition rate by w defined from mi- 
croscopic properties of the system (for instance, according 
to given Hamiltonian of interaction and the Fermi's golden 
rule) . In this case the kinetic equation is of the form 



h:,{t)-w{n^{t)-n:,{t)} = 0, 



(7) 



where the dotted symbol means the first-order time deriva- 
tive. The relaxation function for the two-state system 
reads then 4>D{t) = 1 — 2n_|_(i) = 2n-^{t) — 1 = exp{—2wt). 
It is easy see that the steady state of the system corre- 
sponds to equilibrium with n^(cxD) = nj^(oo) = 1/2. How- 
ever, if the dipoles interact with their environment, and 
the interaction is complex (random), their contribution in 
relaxation already will not result in any exponential decay. 

Assume that the interaction of dipoles with environment 
is taken into account with a help of subordination in time. 
Take the process 5*4 as a subordinator. It accounts for the 
amount of time when a dipole does not participate in mo- 
tion. The ratio of dipoles in the state t and another in the 
state I is subordinated by the process 5*4. The equation 
describing the two-state system takes the form similar to 
Eq. ^ , but the derivatives of first order become fractional 
of order < a < 1 determined by the index of the pro- 
cess St inverse to the a-stable process T'(r). This leads to 
the Cole-Cole (CC) relaxation (see [i[7] for details). The 
relaxation function for the two-state system satisfies now 
the following equation 



'CC 



it) 



r 



r(i - a) 



-Wp</'CCW: 



with the initial condition 0cc(O) = 1- Here we use the 
Riemann-Liouville definition of fractional derivative |10j . 



p-2 



Subordination model of anomalous diffusion . . . 



namely 



© gives 



5" , , d" 
— x t) = — 



1 



x{t) 



r(n - a) Jo (i - t)"+i-" 



dr 



The constant cj^ characterizes 



where n — 1 < a < n 
the transition rate from microscopic properties of the sys 
tem. The relaxation function reads (pcc{t) = Ea{—ujpt°'), 
where Ea{z) = X]fc°=o 2'^/r(l + ka) is the one-parameter 
Mittag-Lefher function [TTj . 

For the experimental study the frequency-domain rep- 
resentation of the latter function 



(p*{io) oc / e 



dt J 



(8) 



is of interest. It is well known that the complex dielec- 
tric susceptibility = x'{'^) ~ *x"('*^) oc '^*{oj) of most 
dipolar substances demonstrates a peak in the loss com- 
ponent x" {'^) a characteristic frequency ujp. The CC 
susceptibility is 



1 



1 + (iw/wp)" ' 



< a < 1 



(9) 



With reference to the theory of subordination the CC law 
shows that the dipoles tend to equilibrium via motion al- 
ternating with stops so that the temporal intervals be- 
tween them are random. However, there are other well- 
know laws of relaxation, in particular, the Cole-Davidson 
(CD) and Havriliak-Negami (HN) ones. The description 
of the laws requires a modification in the theory of subor- 
dination. 

It should be noticed that the physical mechanism un- 
derlying the nonexponential relaxation can be described 
as a diffusive limit of continuous time random walks. The 
resulting relaxation patterns are connected not only with 
stochastic features of the jumps and the inter-jump times 
themselves, but also with a stochastic dependence between 
them. In the framework of Linear Response Theory the 
temporal decay of a given mode k, representing excitation 
undergoing diffusion in the system under consideration, 
is given by the inverse Fourier transform of the diffusion 
front [12]. 

Subordination by the process T{St). The sub- 
ordinator T{St) results in stretching of the real time t. It 
will underline scaling properties in short and long times, 
respectively. The useful feature is just observed in CD and 
HN relaxation. Consider the process T{St) as a subordi- 
nator to exponentially decreasing states. Let it be indexed 
by 7, i. e. the process is obtained from a 7-stable random 
process. The relaxation function for the two-state system 
takes the form 



r(7)r(i-7) 



dz . 



(10) 



Here the subscript CD is not by chance. It shows a di- 
rect connection of the relaxation function with the Cole- 
Davidson law II . In fact, the one-sided Fourier transform 



XCD(t^) 



[1 + iuj/ojpp ' 



< 7 < 1 



It should be mentioned that the theory of subordination 
suggests also one more scenario leading to the CD relax- 
ation. It is based on the inverse tempered a-stable process 
(see [13] in more details). 

Relaxation from compound subordinators. — 

The CC and CD relaxations are only special cases of 
the more general HN law. To get that law, the opera- 
tional time St = Sa(t) of the CC diffusion mechanism has 
to be modified [8] by means of coupling between jumps 
and interjump times in the underlying continuous time 
random walk scheme. In other words, the temporal de- 
cay of a given mode, representing excitation undergo- 
ing diffusion in the relaxing system, will be characterized 
by short- and long-time power laws with different frac- 
tional exponents (as in the HN case) only if the anoma- 
lous diffusion scenario is based on a compound opera- 
tional time. To construct such an operational time, de- 
note conveniently the processes T~{St) and T{St) as 
and X^^ , respectively. They corresponds to the under- 
and overshooting subordination scenarios [S]. Next, we 
can write Z^,^{t) < Salt) < Z2,^{t) for t > 0, where 
Z^^^it) = x'^[So.it)l ZO^it) = XO[S^{t)]. The over- 
shooting subordinator leads (stretching the operational 
time Sa{t)) to the HN relaxation in the form 



z-1 (z - 1)^^ 

r(7)r(i-7) 



dz. (11) 



By direct calculations of the Fourier transformation ^ 
the susceptibility reads 



Xhn(w) = 



1 



[1 + (iw/a;p)"]T 



< Q!,7 < 1 . 



The above approach demonstrates clearly a success in the 
probabilistic treatment of the observed relaxation laws. 
Therefore, we continue our analysis as applied to the un- 
dershooting (compressing Sa{t)) subordinator X!:^[Sa{t)]. 
In this case we obtain a new relaxation function 



q^{t) = E^ - {ujpt)'^z 
Jo 



z^i-^ (1 - z)-' 

r(7)r(i-7) 



■ dz , 



(12) 



which can be identified as a special case of the generalized 
Mittag-Lefiler function [13] 



a, 13 



(x) 



(7,fc)a;'= 



^ T(ka 

k=0 ^ 



a, /3 > 0, 



where (7, k) = 7(7-1- 1) (7 + 2) . . . (7 + fc — 1) is the Appell's 
symbol with (7, 0) = 1, 7 7^ 0. To avoid any confusion, 
it should be mentioned that the two-parameter Mittag- 
Leffler function Ea jj{z) = ^ z'' /r{ka + P), more common 
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Fig. 2: (Color online) Imaginary term of the dielectric suscep- 
tibility for the Havriliak-Negami, Cole-Cole relaxation and the 
function x'i^) (o^ = 0.9, 7 = 0.8). 

in literature pT, is a special case of the function E^, ^(2;) 
with 7 = 1. From the series expansion of the ordinary 
Mittag-Lefller function Ea{z) it is easy to check by di- 
rect calculations of Eq.(ini) that (f){t) = El j^{-{u!pt)"-). 
This type of the relaxation function has been derived in 
the continuous time random walk framework by Jurlewicz 
and Weron in |8]. Using relation with the Mittag-Lefher 
function we may write now the kinetic equation for (jl2p 
in the pseudodifferential equation form 



where d^/dt" is the Riemann-Liouville fractional deriva- 
tive [To], and (j){0) = 1 the initial condition. Taking 
the Fourier transform ([5]), we get the susceptibility corre- 
sponding to ([T^ in the form useful for fitting the dielectric 
spectroscopy data 



[1 + {tuj/ujp)-°'p 



0<a,7<l. (13) 



This result points to the following relationship with the 
HN relaxation function (see Fig. [5]) 

X*{uj) 1 - (itj/a;p)"''xHN(w) . 

For 7 — > 1 the probability density of tends to the 
Dirac J-function. It is easy to check lim^^i ip{t) = (pccit), 
and the kinetic equation for (j){t) takes the Cole-Cole form 
mentioned above. 

Interpretation of experimental data. It is in- 
teresting to compare xhn('^) and X*i'^) in tti^ framework 
of experimental data results. The fractional two-power 
relaxation dependencies 

X(i^) ~ {iuj/ujpY''^^ for > ojp , 

^x(w) = x(0) - xi^) ~ (ii^/ujp)'" for cj < 



Fig. 3: (Color online) Log-log scale fitting example of 
the dielectric spectroscopy data obtained for gallium doped 
Cdo.99Mno.01Te mixed crystal (see more details in the text). 
The straight and dotted lines represent fits to susceptibility 
data points for T = 77 K. The parameter values a ~ 0.97, 
7 f» 0.78, x(0) ~ 1.0165, ujp/2n ^ 507 kHz are estimated by 
the method of least squares. In the frequency range 10*- lO'^ Hz 
one observes superposition [5] of direct current and loss pro- 
cesses indicated by a deviation of x"(.^) data points from the 
theoretical curve. 



are widely observed for many materials in experimental re- 
laxation studies [HE]- The exponents n and m fall in the 
range (0, 1), and ujp denotes the loss peak frequency. The 
HN relaxation is characterized by the exponents m = a 
and m > I — n — 07, and it fits the so-called typical 
relaxation processes. The function x*{^) demonstrates 
up-down with m = 07 and m < 1 — n = a, and it fits 
the less typical two-power-law relaxation pattern which, 
as shown by the experimental evidence, cannot be ne- 
glected (see e. g. [T1[51[TS] and references therein). Such 
a less typical behavior has been also observed by us in 
gallium (Ga)-doped Cdo.99Mno.01Te mixed crystals [16] 
(see Fig. [S] where sample frequency-domain data mea- 
sured for Cdo.99Mno.oiTe:Ga at 77K is fitted with the 
function p^ ). This material belongs to semiconductor 
of group II- VI possessing deep metastable recombination 
centers. Formation of such centers in Cdi_xMnxTe:Ga re- 
sults from the bistability of Ga dopant which makes this 
mixed crystal as an attractive material for holography and 
high-density data storage (optical memories). 

In the language of subordinators this means that the 
process ^^(t) makes a rescaling for small times, and the 
process X^^(t) turns on a similar rescaling for long times. 
As for < a, 7 < 1, in the case of HN relaxation the dec- 
lination of the imaginary susceptibility x"{'-^) foi' fre- 
quencies will be greater than for high frequencies, whereas 
the less typical relaxation shows an opposite relation. 

The original HN relaxation [Tl[2] with exponents < 
a,7 < 1 satisfies m > 1 — n. Its modified version [15], 
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proposed to fit relaxation data with power-law exponents 
satisfying m < 1 — n, assumes < a, 07 < 1. Unfor- 
tunately, the HN function with 7 > 1 cannot be derived 
within the framework of diffusive relaxation mechanisms. 
Only for 7 < 1 the origins of the HN function can be 
found within the fractional Fokker-Planck [T^ and con- 
tinuous time random walk 8 approaches. The approach 
considered above includes all the data in one, mathemat- 
ically unified approach. 

Conclusions. — We have discovered a novel law 
of relaxation accompanied by the well-known Havriliak- 
Negami function. Earlier the development of the theory 
of nonexponential relaxation went behind the fitting of ex- 
perimental data. Now the subordination approach allows 
one to explain not only the well-known experimental laws 
of relaxation, but it makes a prediction of other adequate 
models useful for experimentalists. 
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